Crossover dynamics of dispersive shocks in Bose-Einstein condensates characterized 

by two and three-body interactions 
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We show that the perturbative nonlinearity associated with three- atom interactions, competing 
with standard two-body repulsive interactions, can change dramatically the evolution of ID dis- 
persive shock waves in a Bose-Einstein condensate. In particular, we prove the existence of a rich 
crossover dynamics, ranging from the formation of multiple shocks regularized by coexisting trains 
of dark and antidark matter waves, to ID soliton collapse. For a given scattering length, all these 
different regimes can be accessed by varying the number of atoms in the condensate. 
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Introduction. — Bose-Einstein Condensation (BEC) 
has been successfully described in the mean-field limit 
by the Gross-Pitaevskii equation (GPE), which embodies 
quantum two-body interactions (s- wave scattering) into 
its cubic nonlinear term [l|. Though not compulsory for 
(ideal, i.e. interaction-less) condensation, such term sus- 
tains several, remakable and otherwise unobservable, co- 
herent phenomena such as soliton formation |2[L quantum 
shocks [3, |i| , collapse [5j , nonclassical states [6[ , time re- 
versal, quantum chaos and Anderson localization |7]. The 
exceptionally huge tunability (in magnitude and sign) of 
the scattering length a through Feshbach resonances 
have made the area of ultracold atoms a prolific ground to 
investigate new quantum and nonlinear physics. Partic- 
ular interest, in this context, was stirred by higher-order 
few-body interactions [9[ boosted also by the observation 
of resonant Efimov states [10]. Most efforts in this area 
are presently aimed at assessing the impact of few-body 
recombination loss coefficients which limits the lifetime 



of condensates [ll|. However, almost unexplored is also 
the effect of three-body elastic collisions, which naturally 
arises as a quintic conservative (Hamiltonian) term in the 
higher-order expansion of the GPE H, ElLIiai whose im- 
pact clearly grows as higher densities are being reached 
in experiments [14[. It is therefore of paramount impor- 
tance to assess whether the three-body interactions can 
lead to a clear signature in terms of qualitative new sce- 
narios, especially if the latter becomes accessible in the 
perturbative limit. 

In this letter we address this question with reference to 
the dynamics of dispersive shock waves (DSWs), which 
form in a repulsive BEC when the kinetic spreading 
regularizes the tendency driven by the nonlinearity to 
form multivalued wave fronts. At variance with earlier 
0,11 EH and more recent studies on DSWs which 
are all based on the standard GPE, we investigate the 
case characterized by an attractive three-body nonlin- 
ear correction to the repulsive s-wave scattering. This 



regime is expected to be achievable by exploiting res- 
onance tuning for bosons and turns out to be relevant 
also for superfl uid fermions [20| . Unlike the case of com- 
peting nonlinearities of the same sign (where three-body 
contributions merely strengthen the nonlinearity at high 
densities), we find that the presence of the quintic term of 
different sign leads to a rich phase-diagram of the wave- 
breaking phenomenon, with novel regimes that uniquely 
identify the contribution of the three-body interaction. 
Remarkably, for a given scattering length, the crossover 
between different dynamics can be totally controlled by 
changing the number of atoms in the condensate. Due to 
the ubiquitous nature of DSW, these findings are also of 
fundamental interest in other areas of classical physics, 
including nonlinear optics 21], electronic systems [22j, 
and granular chains [23|. 

Model and general shock scenario. - We start from 
the following dimensionless cubic-quintic GPE, describ- 
ing the free evolution in ID (after trap is released along 
x) of a BEC ruled by two and three-body interactions: 



dip e 2 d 2 ip 
dt + ~2~dx~ 2 



(1) 



having introduced dimensionless spatial x = * yj mNuj / h 
and temporal t = TNuj/e coordinates, with e = 
h/muRTF <C 1 (Rtf is the Thomas- Fermi radius), and a 
rescaled wave function i/j = \I> \/ 'Anha/muuN , being m the 
single boson mass, N = J \^\ 2 dX the number of particles 
and u the transverse BEC trap frequency. The parameter 
a/huj = Ngs/\g2\ 2 (with a > 0) uniquely determines the 
properties of the system in terms of the number of parti- 
cles and the relative strength of the three-body term g% 
Q over the two-body term #2 = ^h 2 a/m. At variance 
with the cubic GPE system (a = 0), which is completely 
characterized by the scattering length a, this normal- 
ization clearly shows that the high-order BEC dynamics 
lives in a two dimensional phase space encompassing both 
a and N. This will permit to explore different scenarios 
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FIG. 1. Color Online. Dispersionless evolution for pj = 1.5: 
(a-b) as obtained from Eqs. (5-7) with a = 0.15 (a), a = 0.3 
(b); (c-d) numerical results from Eqs. (2) comparing the case 
a = 0.4 (c) with the case of a standard GPE (a = 0) with 
positive cubic nonlinearity (d). 



even for a constant scattering length. We study the gen- 
eral shock dynamics in the dispersionless (or hydrody- 
namic) limit of Eq. (pQ), which can be obtained by apply- 
ing the Madelung transformation ijj = v /pexp[^ J dxu) 
for e -> 0: 

dp d{up) du du df 

dt dx dt dx dx 

where f(p) — p — ap 2 , and p, u act as density and veloc- 
ity of the fluid, respectively. The impact of the quin- 
tic term can be understood by studying the decay of 
an initial jump in the density, which can be easily real- 
ized experimentally. This amounts to solve the Riemann 
problem associated with the initial value u(x, 0) = 0, 



p(x,0) 



(p_ — Pq)Q(x) [Q(x) is the Heaviside 



function], pj (pq > p^) being the initial constant den- 
sities across x = 0. We solve Eqs. (j2j) by deriving the 
following Riemann invariants A± = u ± R(p): 



R(p) = " 2ap) + 



cos 1 (1 — 4ap) 



2V2a 



(3) 



dt 



which transform Eqs. ([2]) into the diagonal form 

V±^r^r = 0, where V± = u ± \J p(l — 2ap) are the Rie- 
mann characteristic eigenvelocities. As long as V± are 
strictly real, i.e. 



a< l/(2p ), 



(4) 



the evolution of X± can be given in terms of the self- 
similar variable £ = x/t through the following simple 
wave solutions of the system (characterized by one of the 
two Riemann variables remaining constant) 



R{p+) - R(p) - y/p(l - 2ap) = C, 

R(p) - R(po) + ^p{l - 2ap) = C (5) 
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FIG. 2. Color Online. Global phase diagram of the shock 
dynamics in the parameter plane (po = Po , ot) [pj = 1] ■ 

The two simple waves ([5]) are connected by an interme- 
diate region with constant values p = pi and u = u\. 
The latter can be found by the matching equations 
m = [R(pt) ~ R(Po)}/2 and 2R( Pi ) = R(p+) + R(p^). 
The borders Xi± of the intermediate region are then cal- 
culated by substituting the density pi into Eqs. (j5j): 



^ = Ci± = R(Pi) - R(Po) ± y/piO- ~ 2 <* Pi ). (6) 

Conversely, the simple wave external edges x e ± = ( e ±t 
are calculated by matching Eqs. ((5]) with the asymptotic 
(input) values p^ 1 , thus obtaining 



p+(l-2ap+), 



t 



Vl - 2a. (7) 



Equations dS])-© predict a rich scenario for the shock 
dynamics, as illustrated in Fig. QJi-b. Henceforth, with- 
out loss of generality, we take pj = 1 and p = Pq • For 
small values of a [see Fig. [TJi], wave-breaking occurs only 
for x > 0, where the velocity £i+ of the point x i+ is larger 
than £e+, leading to a multivalued region, whose dis- 
persive regularization leads to a DSW. However, a more 
complex dynamics is observed when a is increased [see 
Fig. [TJd]: in this case not only Xi+ but also x e - moves 
sufficiently fast and induces wave breaking for x < 0, 
thus producing a double shock in the dynamics. The 
crossover between these two regimes is given by the con- 
dition C e _(a*) = Q_(a*). It is worth noting that a 
double shock have also been predicted in a BEC flow- 
ing through a penetrable barrier [17|, though both the 
underlying mechanism and the regularization are com- 
pletely different in the latter case. A third regime settles 
in when Eq. (|4|) is violated, resulting in imaginary eigen- 
velocities V± [Eqs. ([2j) are no longer hyperbolic]. In this 
case, solutions ©-(O are meaningless and we resort to 
numerical integration of Eqs. (j2j). Our analysis show 
(Fig. [it) that the system evolution is characterized by 
the generation of two opposite velocities (Fig. [Th, u in 
the inset), which compress the wavefront and generate 
a cusp-like singularity at x ~ 0~. This behavior orig- 
inates from a three-body dominance over the two-body 
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FIG. 3. Color Online, (a) Velocity of DSW edges 0,0, and 
rarefaction edges f e _, 0- vs - a 'i ( D-c ) Color level plot of 
density evolution p(x,t) from Eq. (1) for e = 0.01 and (b) 
a = 0, (c) o; = 0.1; For comparison the edges of the DSW 
xz, xt from Eq. ©, x e _, x%- of the rarefaction wave from 
Eqs. ©-(0) are reported, (d) Snapshots p(x,t = 8) from 
(b-c) showing the DSW dark soliton train for a = 0, 0.1. 



term, as demonstrated by comparing Fig. []J to Fig. [T]i, 
which shows the evolution of the same input launched in 
a standard GPE with positive nonlinearity [Eq. (pQ) with 
a = and x —> ix]. The two wave-breaking are qualita- 
tively identical, which can be understood by considering 
that the GPE with attractive two-body nonlinearity ex- 
hibits complex eigenvelocities, in the same way as Eq. 
(P) for a > 1/2 po. Despite such similarity, the long-term 
dynamics (regularization) in the presence of three-body 
nonlinearities will be totally different, as discussed in the 
following paragraphs. Figure [2] summarizes the results 
of the dispersionless analysis in a wave-breaking phase 
diagram. In particular, the two curves a* = a*(po) and 
a = l/(2/?o), divide the phase space into three distinct re- 
gions, each characterized by different wave breaking sce- 
narios deepened below. 

Real eigenvelocities: tuning the shock dynamics. 

We study the shock regularization for a < a* by 
exploiting Whitham theory of modulation, in the form 
suitable for nonintegrable systems [24[. We summarize 
only the outcome of this approach, while deferring the 
(involved) mathematical details to a successive paper. 
When a < a*, the multivalued region for x > (Fig. 
[1]) is regularized by the formation of a single DSW such 
that a modulated cnoidal wave (or dark soliton train) 
appears within a shock fan limited by a leading x\ and 
a trailing x t edge. By matching the dispersionless limit 
with the Whitham equations [24J , we are able to find the 
velocities of the edges = x\/t and ( t = x t /t (replacing 
the hydrodynamic estimates £ e+ , i n the form: 



= (l-2a)(2 7/ -- 
C* = R(pi) - R(l) + jtVpi^-^pi), (8) 
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FIG. 4. Color Online, (a) Snapshots of density p{x) obtained 
from Eq. (1) with a — 0.3, e = 0.01, and a input jump po = 
1.5; (b)-(c) zoom of the boxed evolution in (a), showing the 
opposite type of regularization taking place in the dynamics. 



with ji = 7(1) and j t = j(pi) arising from the solution 
of the differential equation: 



dj _ (l + 7)(l + 27-87a) 
~dp ~ 2p(27 + l)(2ap-l) ' 



(9) 



integrated in p G [l,Pi], with the initial condition 7(1) = 
1 for 7/, and j(pi) = 1 for j t . Figure [3^ displays the 
DSW edge velocities £/, ( t [calculated by integrating Eq. 
(j9j)] and those of the rarefaction wave £ e _, 0- [from Eqs. 
©-(jTj)] as a function of a. Increasing the three-body 
contributions results into a nearly linear decrease of all 
relevant velocities £/, 0, Ce- 7 0-? with a consequent re- 
duction of the extension of both simple waves and the 
shock fan. It is worth emphasizing that the two edges of 
the DSW behave in a markedly different way: the leading 
edge (1 is strongly influenced by a, while the variation of 
is much smaller. As a consequence, even when a is 
increased by a small factor, the shock dynamics is sig- 
nificantly affected in terms of both shock fan and shock 
overall angular direction. To verify these predictions, we 
resort to numerical integration of Eq. (pQ). As shown in 
Fig. 03>c, a remarkable agreement is obtained between 
theory and numerical simulations. The snapshot in Fig. 
(3jd shows that, despite a small variation of a, a substan- 
tial reduction of the shock fan extension occurs (nearly 
a factor of two for Aa = 0.1) accompanied by a reduced 
average velocity (center of mass closer to x = 0). In 
summary, in the regime where a single shock is formed 
for x > 0, all the features of the shock dynamics are sig- 
nificantly varied by slightly perturbing the strength of 
the three-body interaction terms a. 

Real eigenvelocities: multi-shock generation through 
antidark solitons. - When a* < a < l/2po, simulations 
show the occurrence of a second DSW for x < (Fig. 2J, 
in complete agreement with the hydrodynamic analysis. 
In this case, however, the DSW for x < differs sub- 
stantially from the dynamics discussed above for x > 0. 
The latter, in fact, is characterized by a cnoidal wave 
composed by a train of dark-like oscillations (see Fig. [3]), 
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FIG. 5. Color Online. Snapshots of density p vs. x from Eq. 
(1) for a = 0.4 and e = 0.01: (a) shock regularization; (b) 
antidark collapse instability. 



dramatically alter the breaking scenario in repulsive 
BEC. The generation of multiple DSWs involving anti- 
dark soliton trains, as well as ID collapsing dynamics are 
new peculiar behaviors that can be controlled by means 
of the BEC number of atoms. These results are expected 
not only to foster new perspectives in BEC physics, but 
also to stimulate novel experiments in nonlinear optics, 
where the quintic terms account for saturation of clas- 
sical Kerr nonlinear it ies. We acknowledge funding from 
PRIN 2009 project (No. 2009P3K72Z). 



owing to the condition p(x t ) > p(xi) (see also Fig. QJi). 
The DSW, more specifically, originates from a series of 
oscillations that start from the leading edge x\ (where the 
shock matches the linear wave background at p = 1) and 
culminate into a dark soliton at the trailing edge [Fig. 
Hfc]. On the contrary, the DSW for x < moves in the 
opposite (backward) direction and has a larger density 
on the leading edge, which yields p(xn) > p(x t i), being 
x\\ and x t i the leading and trailing edge of the DSW at 
x < 0, respectively. As a consequence, the modulated 
wave train that regularizes the shock needs to be com- 
posed by bright entities, as shown in Fig. HJd. In this 
case, the trailing edge is an antidark (bright on pedestal) 
soliton solution of Eq. (1), which has been thoroughly 
investigated recently [25]. Quite remarkably, although 
these solitons are mostly unstable [25j, the average ve- 
locity of the DSW for x < is such that the soliton (and 
hence the DSW) is totally stable on propagation. 

Complex eigenvelocities: ID BEC collapse. — When 
a > l/2po, the wave-breaking observed in the simula- 
tions of Eq. (1) reflects indeed the dominant charac- 
ter of three-body interactions (Fig. EK-b). In the early 
stage, a cusp-like singularity in x ~ _ is generated (Fig. 
[5^), as predicted by Eqs. (j2j). Contrary to the the case 
of the standard attractive GPE, which shows a similar 
early stage behavior [see Fig. [TJi], the three-body dom- 
inance in Eq. (pQ) leads to a completely different long 
term dynamics. In the quintic case, in fact, the singu- 
larity tends to evolve into antidark solitons, due to the 
condition p(xi) > p(x t ). Contrary to the previous case, 
however, here the solitons velocity is nearly vanishing. 
Zero- velocity antidark solitons of Eqs. (1) are always un- 
stable [HI, an d they lead the BEC towards irreversible 
collapse near x = (Fig. [5Jd). We emphasize that such 
a collapsing dynamics, which we have verified to occur 
also for different inputs (e.g. Gaussian on pedestal), is 
a unique features of Eq. (pp). Collapse, in fact, cannot 
be observed in the standard GPE (a = 0) that, owing to 
its integrable nature, does not possess unstable solitons. 
It is also worthwhile remarking the perturbative value of 
the quintic term (a > 0.34) which induces such a dra- 
matic and abrupt change in the dynamics. 
In conclusion, we have shown that the perturbative non- 
linearity arising from three-body elastic collisions can 
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